Breakdown of Stokes-Einstein relation in the supercooled liquid state of phase change 

materials 
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The application of amorphous chalcogenide alloys as data-storage media relies on their ability to 
undergo an extremely fast (10-100 ns) crystallization once heated at sufficiently high temperature. 
However, the peculiar features that make these materials so attractive for memory devices still 
lack a comprehensive microscopic understanding. By means of large scale molecular dynamics 
simulations, we demonstrate that the supercooled liquid of the prototypical compound GeTe shows 
a very high atomic mobility (D ~ 10 -6 cm 2 /s) down to temperatures close to the glass transition 
temperatures. This behavior leads to a breakdown of the Stokes-Einstein relation between the self- 
diffusion coefficient and the viscosity in the supercooled liquid. The results suggest that the fragility 
of the supercooled liquid is the key to understand the fast crystallization process in this class of 
materials. 



I. INTRODUCTION 

Phase-change materials based on chalcogenide alloys 
are attracting a lot of interest due to their ability to un- 
dergo reversible and fast transitions between the amor- 
phous and crystalline phases upon hcatingir— . This prop- 
erty is exploited in rcwriteable optical media (DVD, 
Blu-Ray Discs) and phase change non volatile memo- 
ries (PCM). The strong optical and electronic contrast 
between the crystal and the amorphous allows discrim- 
inating between the two phases that correspond to the 
two states of the memory. The PCM devices, first pro- 
posed by Ovshinsky in the late 1960's£, offer extremely 
fast programming, extended cycling endurance, good re- 
liability and inexpensive, easy integration. A PCM is 
essentially a resistor of a thin film of a chalcogenide alloy 
(typically Gc2Sb2Tc5, GST) with a low field resistance 
that changes by several orders of magnitude across the 
phase change. In memory operations, cell read out is per- 
formed at low bias. Programming the memory requires 
instead a relatively large current to heat up the chalco- 
genide and induce the phase change, either the melting 
of the crystal and subsequent amorphization (reset) or 
the recrystallization of the amorphous (set). 

The key property that makes these materials suitable 
for applications in PCM is the high speed of the trans- 
formation which leads to full crystallization on the time 
scale of 10-100 ns upon Joule heating. What makes some 
chalcogenide alloys so special in this respect and so dif- 
ferent from most amorphous semiconductors is, however, 
still a matter of debate. 

In this paper, we demonstrate by means of atomistic 
molecular dynamics (MD) simulations that the phase 
change compound GeTe displays a very high mobility in 
the supercooled liquid phase down to very low temper- 
atures T leading to a breakdown of the Stokes-Einstein 
relation D oc — (SER) between diffusivity D and viscos- 
ity 77. The breakdown of SER in the supercooled liquid 
phase, which is typical of a fragile liquid^, is actually 



the key to understand the fast crystallization behavior of 
these materials as discussed below. 

Supercooled liquids are classified as fragile or strong 
on the basis of the temperature dependence of their 
viscosity^. An ideal strong liquid shows an Arrhenius 
behavior of the viscosity rj from the melting tempera- 
ture T m down to the glass transition temperature T g . 
On the contrary, in a fragile liquid r\ follows an Arrhe- 
nius behavior only above a cross-over temperature T* 
below which a Vogcl-Tammann-Fulcher (VTF) function 
77 = rj exp( kB {J?-T ) ) * s cus tomarily used to reproduce 
the data with r] Q , E and T D as fitting parameters^. The 
viscosity in a fragile liquid shows a steep rise by approach- 
ing T g which according to the SER would lead to a strong 
decrease in the atomic mobility. On the other hand, the 
self diffusion coefficient D controls both the speed of crys- 
tal growth u and the steady state nucleation rate I ss . 
In fact, classical nucleation theory^ predicts that u oc 
D(l - cxp(-AG/(fc B T)) and I ss oc D exp(-G c /(fc s T)) 
where AG is the free energy difference between the liq- 
uid (or amorphous) and the crystalline phases and G c 
is the formation free energy of the critical nucleus given 
in turn by G c = 167rcr 3 /(3AG 2 ) with a interface energy. 
AG is the driving force for crystallization that decreases 
by increasing temperature and finally vanishes at T m . In 
phase change materials, due to the breakdown of SER, 
the diffusivity can be very high just above T g in spite of a 
large viscosity. Consequently D can reach high values at 
temperatures much lower than T m where a large driving 
force for crystallization is present. 

A breakdown of SER in GST has actually been sug- 
gested by recent measurements of the crystallization rate 
by Orava et al.— . Still, further evidences of the break- 
down of SER are needed to support the crucial assump- 
tions made by Orava et al. to infer the viscosity from 
their only data of differential scanning calorimetry. 

To this aim we performed MD simulations by using a 
classical interatomic potential^ we generated by fitting a 
large database of density functional energies by means of 
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the Neural Network (NN) method introduced in RefAi. 
So far we have restricted ourselves to the binary com- 
pound GcTe which is also under scrutiny for applications 
in PCM and which shares most of the properties with the 
more commonly used GST— i 12 i 13 . The NN potential dis- 
plays an accuracy close to that of the underlying density 
functional theory (DFT) framework at a much reduced 
computational load that scales linearly with the size of 
the system. It allows us to simulate several thousands 
of atoms for tens of ns, which is well beyond present-day 
capabilities of DFT molecular dynamics. 



II. COMPUTATIONAL DETAILS 

The NN interatomic potential of GeTe was obtained 
in Refi^ by fitting a huge database of DFT energies 
by means of the method introduced by Behler and 
Parrinello^i. The database consists of the total energies 
of about 30000 configurations of 64-, 96-, and 216-atom 
supcrcells computed by employing the the Pcrdcw-Burkc- 
Ernzerhof (PBE) exchange and correlation functional^ 
and norm conserving pseudopotentials. The NN poten- 
tial displays an accuracy close to that of the underly- 
ing DFT-PBE framework whose reliability in describing 
structural and dynamical properties of GeTe and other 
phase change materials has been validated in several pre- 
vious works by ou r 10 i 13 i 15 and other groups^. 

The simulations were performed with the NN code 
RuNNcr— and a 4096-atom cubic supercell by using the 
DL_POLY— code as MD driver. The time step was set to 
2 fs, and constant temperature was enforced by a stochas- 
tic thermostat^. 

It turned out that in order to reproduce the equilib- 
rium density of the liquid at T m , an empirical van der 
Waals (vdW) correction had to be added to the NN po- 
tential. This was done by using the scheme proposed by 
GrimmoiS with the sq parameter tuned to 0.55 to repro- 
duce the experimental equilibrium volume of the liquid 
at T„^2. The experimental equilibrium volume of the 
amorphous and crystalline phases are instead well repro- 
duced by the NN potential without the need of the vdW 
interaction-^. The inability of the NN potential in repro- 
ducing the equilibrium volume of the liquid can be traced 
back to the presence of nanovoids in the melt^. In the 
liquid the nanovoids can coalesce and increase in size by 
decreasing the density which results into a reduced tensile 
stress upon expansion. This effect is hindered by vdW 
interactions. Nanovoids are also present in the amor- 
phous phased, but their distribution can not change with 
temperature because of the low atomic mobility in the 
amorphous phase. The calculated linear thermal expan- 
sion coefficient (a=jyj^) of the liquid at T rn with the 
NN+vdW potential turned out to be 4.73 • 10~ 5 K" 1 to 
be compared with the experimental value of 3.73 • 10~ 5 

K -120 

The added vdW interaction acts just as a volume de- 
pendent term in the equation of state of the liquid but 



it is not included in the MD simulations discussed be- 
low. DFT calculations of the self-diffusion coefficient at 
two selected temperatures were performed by molecular 
dynamics simulations within the Born-Oppenheimer ap- 
proach by using the code CP2K— ; Kohn-Sham orbitals 
were expanded in a Triple- Zeta- Valence plus Polarization 
Gaussian-type basis set and the charge density was ex- 
panded in a plancwave basis set with a cut-off of 100 Ry 
to efficiently solve the Poisson equation within the Quick- 
step scheme^. Goedecker-type pseudopotentials^ with 
four and six valence electrons were used for Ge and Te, 
respectively. Brillouin Zone integration was restricted to 
the r point of the 216-atom supercell. The same frame- 
work was used in our previous DFT molecular dynamics 
simulations of GeTe^^. 



III. RESULTS 

To study the properties of the supercooled liquid, we 
first assessed the ability of the NN potential, and thus 
of the underlying DFT-PBE framework, to reproduce 
T m . The melting temperature was computed by means 
of thermodynamic integration^ that yielded T m =1001 
K very close to the experimental value at normal pres- 
sure of 998 K2£. To obtain T m , we first computed the 
difference in the Helmholtz free energy F between the 
NN system and a reference system for which an analytic 
expression for F is known, at a given temperature T" and 
density p' . Namely 

F NN (T',p')-F ref (T',p')= f dX(U(X)), (1) 

Jo 

where the average is taken over a MD simulation with 
the potential U(X) = XUnn — (X — l)U re f. The temper- 
ature and density were set to the experimental values at 
the melting point at normal conditions^. The reference 
system was chosen as an Einstein crystal for the solid and 
a Lennard- Jones fluid^ for the liquid. In the next step, 
the chemical potentials were evaluated by integrating the 
free energy as a function of density starting from p— 

Vnn(T>,p) = ±F NN {T>,p>)+b{T'){lnjr + 1)+ 

+ + c(V)(2p p>) U 

where N is the number of particles in the simulation 
cell, and the parameters a(T'), b(T') and c(T') were ob- 
tained by fitting the pressure dependence of the density 
using 

P(T',p)=a(T') + b(T')p + c(T')p 2 . (3) 

By equating the chemical potential of the two phases one 
obtains a transition pressure of -0.44 GPa at the chosen 
temperature T'=998 K. From the calculated Clausius- 
Clapcyron equation (dT/dP=6.85 K/GPa from the cal- 
culated AS=AE /T and AV on the theoretical melting 
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line at T=998 K) we then obtained the theoretical melt- 
ing temperature at normal pressure which is T m =1001 
K. 

We then analyzed the properties of the supercooled liq- 
uid below T m by computing independently 77 and D in 
microcanonical MD simulations. The volume of the su- 
percooled liquid was scaled with temperature according 
to the thermal expansion coefficient (a is little dependent 
on temperature) given in Sec. II. We scaled the temper- 
ature from 1000 K to 500 K in eleven steps. At each step 
the system is equilibrated for 25 ps with the thermostat. 
Overall the system is thus quenched from 1000 K to 500 
K in 250 ps. At each temperature statistical averages 
are then collected on microcanonical simulations up to 2 
ns long for the calculation of the viscosity as discussed 
below. 

We first computed D from the atomic mean square 
displacement on the time scale of 50 ps on which the 
system does not crystallize at any temperature we con- 
sidered. The values of D as a function of temperature are 
reported in Fig. [T^,. D is still of the order of 10~ 6 cm 2 /s 
at the lowest temperature of 505 K considered here; it 
follows an Arrhenius behavior from T m to 505 K. The 
activation energy is 0.220 ± 0.002 eV, a value much lower 
than the activation energy of 1.76 eV obtained from the 
Arrhenius dependence of viscosity in GST measured in 
the temperature range 333-373 K (probably below T 9 )2i. 
The ratio between the self-diffusion coefficient of Ge and 
Te {D(j e / Dxe) increases by decreasing temperature as 
shown in Fig. [5J The values of D obtained from the NN 
simulations were also validated by direct DFT molecular 
dynamics simulations at few selected temperatures with 
a small 216-atom cell (cf. Fig. Q}i). The DFT result at 
1000 K is equal to that previously obtained in Ref^ with 
the same cell and functional used here. We checked that 
the change of volume with temperature has a little effect 
on the diffusion coefficient as shown in Fig. QJ1 which also 
reports the values of I? as a function of temperature once 
the density is fixed to the value at the melting point. 

We then computed r\ between T m and a temperature 
T*=700 K which turned out to be our crossover temper- 
ature by means of the Green-Kubo (GK) formula^ 

v= k^fJ K^KfC ))^ ( 4 ) 

where a xy is the off diagonal component of the stress 
tensor and V the supercell volume. The integral in Eq. 3] 
is converged by restricting the integration time to 60 ps 
due to the decay of the self-correlation function above T* . 
However, long simulation times up to 2 ns are needed to 
converge the average (< .. >) over different initial times 
t = 0. 

Above T* the viscosity can be described by a simple 
Arrhenius (Fig. [TJ)) function with an activation energy 
of 0.17 ± 0.035 eV, very close to the value of 0.2 eV mea- 
sured experimentally for the Geo.15Teo.85 eutectic alloy 
above T^. For the GeTc composition, experimental 



values of r\ are available only at 1000 K yielding ?/=2.59 
mPa-s which is twice as large as our result (cf. Fig. [TJd). 
This discrepancy is not due to the NN potential but pos- 
sibly to limitations of the underlying DFT framework. 
Previous works on GeSe2 have indeed shown that differ- 
ent choices of the exchange and correlation functional af- 
fect the dynamical properties of the liquid phased. The 
viscosity can be computed from the GK formula only 
above T* since at lower temperatures the system crys- 
tallizes spontaneously on the time scale of few hundreds 
of ps which is not long enough to get the value of 77 con- 
verged from Eq. @] In the supercooled liquid, r\ can not 
be defined on a time scale longer than the crystalliza- 
tion time which in GcTc is very short in the temperature 
range 500-700 K. 

We thus attempted to extrapolate rj below T* by a 
VTF-like function with the constraint of matching the 
typical value of 10 15 mPa- s expected at T^i. Unfortu- 
nately, a reliable value of T s for is not available from 
experiments because of the fast crystallization of GcTc, 
and its theoretical estimate from simulations is difficult. 
Experimental data on T g are available for the better glass 
formers Ge x Te^i_ x -j alloys with x=0. 15-0.23^. By a lin- 
ear extrapolation with x of these latter data on T g one ob- 
tains T 9 =511 K for £=0.5 which is probably too high. On 
the other hand, T g is customarily assumed to be slightly 
below the crystallization temperature, that is about 450 
K in GeTe^L. We then used the function proposed in22 
that allows fitting r\ over a wider range of temperatures 

iogio V(T) = log 10 Vo + (15 - log 10 rj )- 

T ex P [^15-log 10 , )o ) \ T )_ 

where to and rj are fitting parameters. In Eq. [5] 
? ? =10 15 mPa- s at T g . 

The parameter to is is the fragility index of the super- 
cooled liquid defined by the logarithmic derivative of 77 
at T 9 , to = d(logiorj) / d(T g /T) \t=t s - We performed two 
fittings for two values of T g as shown in Fig. [3] the first 
with T s =450 K that yields m=lll (logior] =-0.18), and 
a second one with a somehow lower temperature T s =400 
K which yields a very similar value of to=104 (logior] =- 
0.15). Similar results are obtained by using the modified 
VFT function proposed in Ref<22. Our data on viscosity 
above T* are thus consistent with a high fragility of the 
supercooled liquid. For sake of comparison we remark 
that m=20 in silica which is a typical strong liquid while 
to=191 in PVC which is a typical fragile liquid^. Un- 
fortunately, due to the lack of data on 77 below T* and 
the uncertainties in the value of T ff , we can not assign 
accurately the degree of fragility. Nevertheless, even for 
a value of to as large as 100, the viscosity rises too steeply 
in the range 500-600 K to be consistent with the calcu- 
lated values of D and the application of SER. In fact, the 
SER actually breaks down at T < T* as discussed below. 

In the hydrodynamic regime when the SER holds it is 
actually possible to estimate the viscosity on the shorter 
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FIG. 1. a) Self-diffusion coefficient D as a function of temperature in the supercooled liquid GeTe calculated from the mean 
square displacement (red triangles). The density is scaled with temperature according to the calculated thermal expansion 
coefficient. The open squares are the results of DFT simulations of a 216-atom cell. The crosses correspond to the values of 
D computed by holding the density fixed to the value at the melting point, b) Calculated viscosity r\ (Green-Kubo formula) 
as a function of temperature in the supercooled liquid. The density changes with temperature. The lines are Arrhenius fits of 
the data that give an activation energy of 0.220 ± 0.002 eV for D and 0.17 ± 0.035 eV for r\. T m is the theoretical melting 
temperature (see text). 
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FIG. 2. Ratio of the self-diffusion coefficients of the two 
species Dgh/Dte as a function of temperature in supercooled 
liquid GeTe. 



time scale of 50 ps by a finite size scaling analysis of the 
self-diffusion coefficient. The viscosity can be obtained 
from the scaling of D with the edge L of the cubic simu- 
lation cell as2£ 



magnitude larger than those obtained from D and the ap- 
plication of the SER (rj = §®^ D , where R is the average 
van der Waals radius of the two species) as shown in Fig. 
|U This inconsistency demonstrates that the SER indeed 
breaks down. We remark that the numerical values of rj 
reported in Fig. @] below T* are not reliable since they 
are obtained from Eq. [6] which is not applicable when 
the SER breaks down. 
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We considered three cubic models with 512, 1728 or 
4096 atoms. By applying Eq. [6] above T*, we obtained 
values for r\ very close to the GK data (Fig. @| and con- 
sistent with the SER. However, when Eq. [5] is applied 
below T* , one obtains values of r\ that are three orders of 



FIG. 3. Viscosity as a function of temperature in supercooled 
liquid GeTe. The lines are fittings with Eq. [5] of the Green- 
Kubo data (triangles, cf. Fig. [lb) above T* under the con- 
straint of reproducing the value of 7?=10 12 Pa-s at T g with 
T 9 =450 K (continuous line) or T 9 =400 K (dashed line, see 
text). A magnification of the data for T>T* is shown in the 
inset. 
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FIG. 4. Viscosity computed from the Green-Kubo (GK) for- 
mula (red triangles), from the scaling of the diffusion coeffi- 
cient with the simulation cell size (stars, Eq. [6}, and from the 
self diffusion coefficient of the 4096-atom cell and the appli- 
cation of the SER (open circles). 

IV. CONCLUSIONS 

We have demonstrated by means of MD simulations 
that the supercooled liquid of the prototypical phase 
change compound GeTe shows a high atomic mobility {D 
~ 1CT 6 cm 2 /s) down to temperatures very close to the 
glass transition temperature. Our calculated values of 
the viscosity as a function of temperature are consistent 
with a high fragility of the supercooled liquid. However, a 
compelling assessment of the degree of fragility would re- 
quire a reliable estimate of T g which is unfortunately un- 
known. The comparison between the calculated self dif- 
fusion coefficient and the viscosity demonstrates a break- 
down of the SER below a crossover temperature of 700 
K. These results support the experimental evidence of a 
breakdown of SER in the similar compound GST inferred 
by Orava et al. from ultrafast DSC measurements^. This 
feature is the key to understand the origin of the high 
crystallization rate in phase change memories. In the set 
process of PCM, the amorphous phase is heated above 
T g . Preliminary simulation results actually show that 
the glass transition experiences a hysteresis, the rapidly 



overheated amorphous phase differing from the super- 
cooled liquid in a narrow range of temperature above T g . 
The overheated amorphous phase displays a breakdown 
of SER as well, although somehow less pronounced than 
in the supercooled liquid. This issue will be addressed in 
a future publication. However, in the set process of PCM, 
the temperature can rise significantly above T g depend- 
ing on the details of the programming pulses. Indeed 
optimized set pulses with peak intensities equal to the 
reset pulse (leading to the melting of the crystal) are im- 
plemented in current PCM devices^. Under these con- 
ditions the crystallization occurs from the supercooled 
liquid phase. In spite of a large viscosity, a high diffu- 
sivity is still possible in the supercooled liquid down to 
temperatures very close to T g because of the breakdown 
of SER. This allows for the coexistence of a high diffusiv- 
ity and a high driving force for crystallization that boost 
the crystallization speed at high supercooling. The crys- 
tallization of the supercooled liquid or of a highly mobile 
overheated amorphous phase might take place in a dif- 
ferent manner with respect to the crystallization of the 
amorphous phase at temperatures below T g which is of 
interest for data retention. The conclusion we can draw is 
that the self-diffusion coefficient in the supercooled liquid 
regime close to T g can not be inferred from the viscosity 
and the application of the SER. Secondly, the measured 
Arrhenius behavior of the diffusivity and/or crystalliza- 
tion speed below T g can not be extrapolated above T g 
in the supercooled liquid. These results are of interest 
also for the refinement of the electrothermal modeling of 
PCM devices^. 
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